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We propose an approximate method for evaluating the importance of non-Born- Oppenheimer 
effects on the quantum dynamics of nuclei. The method uses a generalization of the dephasing 
representation (DR) of quantum fidelity to several diabatic potential energy surfaces and its compu- 
tational cost is the cost of dynamics of a classical phase space distribution. It can be implemented 
easily into any molecular dynamics program and also can utilize on-the-fly ab initio electronic struc- 
ture information. We test the methodology on three model problems introduced by Tully and on 
the photodissociation of Nal. The results show that for dynamics close to the diabatic limit the 
decay of fidelity due to nondiabatic effects is described accurately by the DR. In this regime, unlike 
the mixed quantum-classical methods such as surface hopping or Ehrenfest dynamics, the DR can 
capture more subtle quantum effects than the population transfer between potential energy surfaces. 
Hence we propose using the DR to estimate the dynamical importance of diabatic, spin-orbit, or 
other couplings between potential energy surfaces. The acquired information can help reduce the 
complexity of a studied system without affecting the accuracy of the quantum simulation. 



I. INTRODUCTION 

The nonadiabatic effects play an important role in 
many chemical phenomena and often must be taken into 
account in accurate calculations of molecular properties 
such as spectra or reaction ratesP^l Many nonadiabatic 
quantum (QM) simulations are performed in the diabatic 
(or quasi-diabatic) basis, which offers several computa- 
tional advantages, especially in the vicinity of conical in- 
tersections of adiabatic potential energy surfaces (PESs) . 
Usually, more PESs must be included to describe a sys- 
tem of interest accurately. One may ask: Which diabatic 
surfaces are important? How much is the result affected 
by neglecting the less important surfaces? The method 
we propose below quantifies the importance of couplings 
between PESs and thus can help in answering these ques- 
tions. 

A direct way to determine whether the coupling of di- 
abatic PESs affects the result of a simulation is to com- 
pute the desired quantity by running QM dynamics for 
both the uncoupled and coupled systems, and to compare 
the results. Here, we propose a more general approach, 
which can give the information about the importance of 
couplings for all observables at the same time. It consists 
in comparing the wave functions by computing the QM 
fidelity, defined a:P 

F QM (t) = \f(t)\ 2 = \(Mt)\Mt))\ 2 - (i) 

In the general setting, ipo(t) and ip p (t) are wave func- 
tions evolved in the unperturbed and perturbed systems, 
respectively. In our case the "unperturbed" represents 
the uncoupled system, "perturbed" the coupled system, 
and ip denotes the full molecular wave function, i.e., it 
includes both nuclear and electronic degrees of freedom. 
Values of -Fqm^) close to unity imply that the perturba- 
tion is not important: the uncoupled Hamiltonian can be 



used in quantitative simulations. Values significantly be- 
low unity imply that the perturbation is important: the 
affected PESs and couplings should be included in the 
simulation. 

Instead of computing fidelity directly from the def- 
inition JT]), one can use the dephasing representation 
(DR),HHZT w hich is an efficient semiclassical approxima- 
tion of fidelity.^ Recently, the DR was used successfully 
to evaluate the ac curac y of QM dynamics on a single but 
approximate PES P 1111 Below, we generalize this method 
to several surfaces and use it to evaluate how the dynam- 
ics is affected by couplings between the diabatic PESs. 

Unlike the computational cost of a direct QM calcula- 
tion, the cost of the DR does not grow exp onentially with 
the number D of degrees of freedomPED Hence the DR 
can be used for many-dimensional systems inaccessible to 
current methods of QM dynamics. An advantage of the 
DR in comparison with mixed QM-classical methods for 
nonadiabatic dynamics, such as the Ehrenfest dynam- 
ics, various surface hopping methodsp^HUl or methods 
in which the classical limit is obtained by the lineariza- 
tion of the QM propagator,^ is that the DR, being a 
semiclassical method, takes the nuclear coherence effects 
into account at least approximately. Other semiclassical 
approaches to nonadiabatic dynamics exist, in cludin g, 
amongst others, the multiple-spawning method d 17 ! 18 ! or 
methods^ based on the Miller-Meyer-Stock- Thoss clas- 
sical electron model. The advantage of the DR is that, 
unlike the majority of semiclassical approaches, the DR 
does not require the Hessian of the potential energy, 
which is the most expensive element of first-principles 
semiclassical dynamics methods (see, e.g., Ref. |20|) . 

Clearly, the above mentioned advantages do not come 
for free: First, unlike most other approaches noted above, 
the DR is not a general dynamical method. It can only 
describe properties expressible in terms of fidelity. Sec- 
ond, the DR in the diabatic basis is expected to work ac- 
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curately mainly when the coupling is relatively weak or 
the nuclear velocities are large, i.e., when the dynamics is 
close to the diabatic limit. Nevertheless, in many cases 
the DR works very well even for strong perturbations 
that lead to completely different classical phase space 
distributions for uncoupled and coupled Hamiltonians. 

Applications for which the DR method is particularly 
well suited include the chemical reactions which proceed 
near the diabatic limit. The photodissociat ion of bro- 
mopropionyl chloride or bromoacetyl chloride^ 21 * 22 ! serve 
as examples. Another possible application would be to 
quantify the relative importance of Hamiltonian coupling 
terms of various origins, e.g., of the spin-orbit vs. dia- 
batic couplings. For instance, in the C1( 2 P) + H2 reac- 
tion the spin-orbit coupling terms clearly dominate over 
diabatic couplings, which can therefore be neglected in a 
simulationPSl 



II. THEORY 

The starting point is the Hamiltonian describing nu- 
clear motion in a molecule, expressed in the diabatic ba- 
sis. To compute the decay of fidelity due to the coupling 
between n PESs, the Hamiltonian n X n matrix is split 
into the uncoupled (i.e., diagonal) part H dlag and the 
coupling (i.e., off diagonal) part AV, 



H 

jjdiag 



jjdiag + A ^ 

T + V diag and AV = V offdiag . (3) 



In Eq. T is the diagonal nuclear kinetic energy ma- 
trix and V dlag contains the diabatic PESs, which are un- 
coupled in H dlag and coupled in H by the elements of 
V° g . (Bold face denotes n x n matrices, hat "denotes 
operators.) 

To derive the DR, one starts from the expression 
for QM fidelity amplitude, applicable to both pure and 
mixed statespl and generalized to the multi-PES setting, 

(4) 



where p is the density operator of the initial state. Gen- 
eralizing the derivation from Ref. I24| fidelity amplitude 
can be written exactly as 

/ QM (*) = Tr J dxp w {x)-(e+^l h e~^ H / h )^x), (5) 

where p w is the Wigner transform of the initial state p, 

[pw] i Ax)=h- D [ de- 



2/ exp 



and x denotes the point (q,p) in the 2 x D-dimensional 
phase space. Approximating the Wigner transform of 
the product of the time-evolution operators, one arrives 



at the DR expression 



/dr(*) = Tr 



d.T pw(^)-e- ,;AS ^°'*) /ri 



AS (x°,t) = f drAVw [x T (x )] , 
Jo 



(6) 
(7) 



where AS (x°,t) is the action at time t due to Wigner 
representation AVw of AV along the trajectory x T of 
H dlag . Note that AVw, AS, and p w are matrix quan- 
tities. If AV contains only diabatic coupling elements, 
then AV = AV(q) and AV w (x) = AV(g). For initial 
Gaussian wave packets, the Wigner function equals the 
classical phase space density which is strictly a probabil- 
ity distribution. 

Since Tr J dx°p w (x°) = 1, it follows from Eq. ^ 
that /dr can be computed as a Monte Carlo average 
(exp(— iAS(x°, t)/h) p ^ with initial conditions sam- 
pled from the Wigner distribution p w (x°) of the initial 
state. However, as AS in Eq. ^ is much smaller than 
the action in other SC methods, the DR alleviates the 
notorious "sign problem." 

If the initial state lies on a single surface n, i.e., if 
[pw] n j = and [pw] J3 = for j ^ n, then no other ele- 
ments of V = V diag + V offdiag than V nn and V nj ,j ^ n, 
enter the calculation. This means that no information 
about other diagonal elements Vjj is used in the DR cal- 
culation. Thus, except for special cases, F DR is expected 
to approximate Fqm accurately only when the detailed 
structure of the remaining PESs does not significantly 
affect the dynamics on V nn . 



III. RESULTS AND DISCUSSION 

Tally's model problems. First, the method is tested 
on three one-dimensional model potentials proposed by 
Tully^ to cover the most important characteristics of 
nonadiabatic transitions. Diabatic and adiabatic PESs 
as well as the coupling terms V12 for Tully's problems 
A and C are shown in Figs. T][2 
found in Ref. Q31 



further details can be 



In all cases, the initial wave packet was a Gaussian with 
approximately 10 % dispersion in momentum and located 
on the lower (problems A and B) or upper (problem C) 
energy diabatic PES in the asymptotic region without 
diabatic couplings. The equations of motion were inte- 
grated until the wave packet left the interaction region. 

The simple avoided crossing model (problem A in the 
original papei!^|) represents the most often encountered 
situation. As can be seen in Fig. [T] the survival proba- 
bility Pi.qm — Trp xl (i) on the initial PES is nearly equal 
to the QM fidelity, suggesting that the fidelity decay is 
caused almost exclusively by the transition to the second 
diabatic PES. Close to the diabatic limit, Fdr accurately 
reproduces Fqm- In practice, the error of comes 
from the intrinsic error of the approximation and from 
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the statistical error due to the finite number of trajecto- 
ries. In this case, the intrinsic relative error is ~ 0.5 % 
when fidelity decays by 10 % to 0.9 and even less for the 
wave packets with higher initial momentum and hence 
fidelity. The relative statistical error remains below 1% 
even when only a single trajectory is used, since most of 
the fidelity decay is due to the transitions to the second 
surface. For slower wave packets, fidelity decreases, and 
the agreement between Fqm and Fdr stays qualitative, 
with Fdr, always decaying faster than -Fqm- Finally, for 
a fixed initial momentum, decreasing the coupling V12 
results in a slower decay of fidelity but the relative error 
of (Fqm — -Fdr) /(I — Fqm) is approximately independent 
of the coupling. Conversely, if V12 is increased substan- 
tially, fidelity initially decays quickly to zero, which is 
captured well by the DR. Subsequently, fidelity may rise 
again, which is usually not well reflected by the DR. 

Tully's system A 
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FIG. 1: Fidelity in Tully's problem A. (a) The diabatic and 
adiabatic PESs and the diabatic coupling, (b) and (c) : Quan- 
tum fidelity Fqm, the quantum survival probability Pi,qm on 
the initial PES, and Fdr as functions of time for two different 
values of initial kinetic energy To. 

In the dual avoided crossing model (problem B, not 
shown), the DR works the best again for high energy 
wave packets. For those the detailed structure of the 
second PES does not significantly affect the motion on 
the first PES. At low energies, where a significant transfer 
of probability density to the diabatic surface V22 and back 
occurs, the DR fails to reproduce the QM fidelity even 



qualitatively. Nevertheless, in this case both Fdr and 
Fqm initially decay almost to zero (although they might 
rise again later), correctly reflecting that the coupling is 
important and should not be neglected. 

A very interesting situation occurs in the extended cou- 
pling model (problem C), where the two diabatic PESs 
Vn and V22 are almost equal (they differ just by a small 
constant shift) but the adiabatic surfaces are well sepa- 
rated due to the coupling. At very low energies of the 
wave packet [Fig. |j(b)], fidelity F goes to zero despite 
that the survival probability Pj on the upper surface re- 
mains close to unity. Therefore, in contrast to previously 
discussed cases, it is not possible to estimate the extent 
of the nondiabaticity of QM dynamics on the basis of 
electronic transitions only. Unlike survival probability, 

Tully's system C 
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FIG. 2: Fidelity in Tully's problem C. 

fidelity describes the nondiabaticity effects on the coher- 
ent nuclear dynamics correctly. It is remarkable that 
the DR describes F so accurately despite that the DR 
dynamics on the diabatic surface ignores the reflection 
of a large part of the QM wavepacket from the upper 
adiabatic surface. The decrease of fidelity is correctly 
captured by Fdr for all energies of the wavepacket. At 
higher energies [Fig. |2](c)], fidelity converges to zero af- 
ter several oscillations, whereas the survival probability 
converges to 1/2. At very high energies (not shown), 
two wave packets moving on the two PESs interfere for 
a long time, Fdr and Fqm agree and oscillate between 
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zero and unity. Surface hopping and Ehrenfest dynamics 
might predict Pi quite well, but since Vu ~ V22, surface 
hopping would incorrectly predict that F w Pi. While 
Ehrenfest dynamics would predict a decay of F below P, 
it is unlikely that it would capture F quantitatively since 
the QM wave packet splits into a faster and slower com- 
ponents whereas Ehrenfest dynamics uses a single mean- 
field surface. One of the reasons why the DR performs 
so well in the extended coupling model is the similar- 
ity of the diabatic PESs in the coupling region. If the 
surfaces are different, e.g., when V22 = cx, the fidelity 
decay at low energies is still well approximated. How- 
ever, at higher energies, Fdr oscillations slowly dephase 
from Fqm , since the DR neglects effects of V22 ■ 

Photodissociation of Nal. We also applied the method- 
ology to the photodissociation of Nal using a two-surface 
model of Engel and MetkP^S (see Fig. |3j). In the orig- 
inal experiment of Mokhtari et alW^ the molecule was 
excited by the light in the 310-390 nm range, which led 
to only a weakly nonadiabatic motion of the wave packet 
on the excited surface. In this regime (corresponding to 
T ~ 0.07 a. u. at the plateau of the PES), far from the di- 
abatic limit, Fdr does not quantitatively reproduce -Fqm- 
However, when the momentum is approximately 10 times 
higher and the dissociation of Nal almost diabatic, Fdr, 
reproduces Fqm rather well. 
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FIG. 3: Fidelity in the photodissociation of Nal. 
Computational details. QM dynamics calculations 



used the first-order split operator method. Classical tra- 
jectories were computed using the first-order symplectic 
Euler algorithm. The time step varied from 0.02 a.u. to 
1 a.u., depending on the system and initial momentum. 

To guarantee convergence, 16384 classical paths were 
used to produce the DR plots. The statistical error of 
/dr for a two-surface system due to a finite number N 
of paths is given by 



4at = ^(( COs2 



(</>)) - (cos((^)) 2 ), 



(8) 



where ip = ASi2(x°,i) /H is the phase accumulated 
along a trajectory and (cos (ip)) = /dr- As demonstrated 
on Tally's model A, where a single trajectory was suffi- 
cient, a much lower number of trajectories than 16384 is 
needed to obtain an accurate estimate of fidelity in cases 
where fidelity stays close to unity. Equation Q implies 
that for a given value of cr s tat and /dr, the number N 
of trajectories needed is approximately independent of 
dimensionality. 

Conclusions. Presented results demonstrate the util- 
ity of the DR in analyzing the molecular QM dynamics 
involving multiple PESs. On one hand, in the nearly dia- 
batic regime, Fdr accurately approximates Fqm- On the 
other hand, in systems far from the diabatic limit Fdr 
decays quickly to zero and thus detects the importance 
of nondiabatic couplings although it may not reproduce 
Fqm accurately. Hence, the method can be used to es- 
tablish the level of nondiabaticity of QM dynamics, with- 
out the need for a QM dynamics simulation. In fact, we 
propose the condition F 1 (instead of the standard 
requirement of high survival probability) as the rigorous 
definition of the diabatic limit (Fig. [2]). Nevertheless, for 
single avoided crossings, e.g., fidelity could be used to es- 
timate the survival probability and hence the branching 
ratios (Figs, [I] and [3j). 

The DR calculation can be performed easily for all sys- 
tems accessible to classical molecular dynamics and for 
which coupling elements Vij,i ^ j, are available. How- 
ever, it remains to be verified how the method will per- 
form in higher-dimensional systems and in systems with 
more than two important surfaces. The first issue was 
partially addressed in Ref. Hill where the DR was ap- 
plied to the two-dimensional photodissociation of CO2, 
confirming that neither the accuracy nor the number of 
classical trajectories needed are significantly affected by 
increased dimensionality. 

At the moment we are exploring a related and comple- 
mentary problem to that of the present paper, namely 
whether the DR in the adiabatic basis can estimate the 
level of nonadiabaticity of QM dynamics near the adia- 
batic limit. 
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